{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Homework 5\n",
    "\n",
    "**Due: 05/07/2020** (Thursday 7th May at 11:59pm).\n",
    "\n",
    "## Instructions\n",
    "\n",
    "+ Develop the code and generate the figures you need to solve the problems using this notebook.\n",
    "+ For the answers that require a mathematical proof or derivation you can either:\n",
    "    - Type the answer using the built-in latex capabilities. In this case, simply export the notebook as a pdf and upload it on gradescope; or\n",
    "    - you can print the notebook (after you are done with all the code), write your answers by hand, scan, turn your response to a single pdf, and upload on gradescope. \n",
    "\n",
    "\n",
    "\n",
    "**Note**: Please match all the pages corresponding to each of the questions when you submit on gradescope. \n",
    "\n",
    "## Student details\n",
    "\n",
    "+ **First Name:**\n",
    "+ **Last Name:**\n",
    "+ **Email:**\n",
    "\n",
    "## Readings\n",
    "\n",
    "Before attempting the homework, it is probably a good idea to review the slides and/or lecture handouts of lectures 19-24 (Inverse problems and Bayesian inference)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "#from autograd import numpy as anp, elementwise_grad as egrad\n",
    "\n",
    "import theano as th\n",
    "from theano import shared, function, tensor as tt\n",
    "import pymc3 as pm\n",
    "\n",
    "import math\n",
    "import scipy.stats as st\n",
    "import scipy\n",
    "\n",
    "import pandas as pd\n",
    "import io\n",
    "import requests\n",
    "\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set()\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 1  - Catalysis problem "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that we used the problem of calibrating reaction rate coefficients in a catalytic reaction as the running example for demonstrating various approaches to solving inverse problems - beginning with the classical approach where this task is posed as the minimization of a misfit function, to the probabilistic approach where the inverse problem is posed as a Bayesian inference task. In this assignment, we will re-visit the catalysis problem (yet again !), this time solving it with ```pyMC```. Working through this assignment should help you get comfortable with probabilistic programming. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Time</th>\n",
       "      <th>NO3</th>\n",
       "      <th>NO2</th>\n",
       "      <th>N2</th>\n",
       "      <th>NH3</th>\n",
       "      <th>N2O</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>500.00</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>30</td>\n",
       "      <td>250.95</td>\n",
       "      <td>107.32</td>\n",
       "      <td>18.51</td>\n",
       "      <td>3.33</td>\n",
       "      <td>4.98</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>60</td>\n",
       "      <td>123.66</td>\n",
       "      <td>132.33</td>\n",
       "      <td>74.85</td>\n",
       "      <td>7.34</td>\n",
       "      <td>20.14</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>90</td>\n",
       "      <td>84.47</td>\n",
       "      <td>98.81</td>\n",
       "      <td>166.19</td>\n",
       "      <td>13.14</td>\n",
       "      <td>42.10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>120</td>\n",
       "      <td>30.24</td>\n",
       "      <td>38.74</td>\n",
       "      <td>249.78</td>\n",
       "      <td>19.54</td>\n",
       "      <td>55.98</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>150</td>\n",
       "      <td>27.94</td>\n",
       "      <td>10.42</td>\n",
       "      <td>292.32</td>\n",
       "      <td>24.07</td>\n",
       "      <td>60.65</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>180</td>\n",
       "      <td>13.54</td>\n",
       "      <td>6.11</td>\n",
       "      <td>309.50</td>\n",
       "      <td>27.26</td>\n",
       "      <td>62.54</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Time     NO3     NO2      N2    NH3    N2O\n",
       "0     0  500.00    0.00    0.00   0.00   0.00\n",
       "1    30  250.95  107.32   18.51   3.33   4.98\n",
       "2    60  123.66  132.33   74.85   7.34  20.14\n",
       "3    90   84.47   98.81  166.19  13.14  42.10\n",
       "4   120   30.24   38.74  249.78  19.54  55.98\n",
       "5   150   27.94   10.42  292.32  24.07  60.65\n",
       "6   180   13.54    6.11  309.50  27.26  62.54"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# load the catalysis data \n",
    "catalysis_url=\"https://raw.githubusercontent.com/PredictiveScienceLab/uq-course/master/lectures/catalysis.csv\"\n",
    "s=requests.get(catalysis_url).content\n",
    "catalysis_data = pd.read_csv(io.StringIO(s.decode('utf-8')))\n",
    "catalysis_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAHwCAYAAAC7T84CAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeZzVZd3/8dcsDCAzA4qTiKCC5eWChBWZ3pi3ZmQmLbfe2e0SKCJ5C1qkKS7kUokL5pqSAhKK2h23plm2kQu3Wlq5pV75K0RlkZFEZhDZZn5/nAMOszBnZr5nzpwzr+fjwQO4rus7fA4fB8/7fL/f61tUX1+PJEmSJKnjinNdgCRJkiQVCgOWJEmSJCXEgCVJkiRJCTFgSZIkSVJCSnNdQIZ6AiOB5cDmHNciSZIkqXsrAXYFngbWN5zIl4A1Eng810VIkiRJUgOHAosaDuRLwFoO8M47a6mrc1v5QtW/fzmrVtXmugxlkT0ufPa48Nnj7sE+Fz573DHFxUXsuGMfSOeUhvIlYG0GqKurN2AVOPtb+Oxx4bPHhc8edw/2ufDZ40Q0uX3JTS4kSZIkKSEGLEmSJElKiAFLkiRJkhJiwJIkSZKkhBiwJEmSJCkh+bKLoCRJkqTtWLduLbW1q9m8eVOra1euLKaurq4Tqso/JSWllJf3o3fvPu063oAlSZIk5bl169ZSU/MO/fpV0aNHGUVFRdtdX1pazKZNBqzG6uvr2bhxA6tXVwO0K2R5iaAkSZKU52prV9OvXxVlZT1bDVdqWVFREWVlPenXr4ra2tXt+hoGLEmSJCnPbd68iR49ynJdRsHo0aMso0stm2PAkiRJkgqAZ66S05G/SwOWJEmSJCXEgCVJkiRJCcmrXQS/c8sTrHxn3dbfV/Yp47rJo3JYkSRJkpT/vnnjItas3dDqumy+/z7uuDGUlJQwd+499OrVa5u5SZNOZ9CgwZx//sUAvPfeWubPn8fChb9lxYrllJdXMGzYcE444WSGDRu+zbG/+MXPufvueSxfvoyBA3fjv/7rZL7whS9m5TVAhgErhLA/8GIzU4fGGBeFEEYDVwEBeBU4L8b4qwbHfwi4CRgNbADmABfGGNt351haJv8RSJIkSdq+TN9XZ/v999KlbzJz5s2cffa3W65hzbuceeYENm7cxIQJ32C//YaxevU7PPDAfZx55gS+850LtwaoRx75PTNmTOfccy9gxIiP8ec/P81VV32fvn37MmrUYVl5DZmewRoGvA0c0Gh8VQhhP+AB4HJgAXAicH8I4WMxxr+l1y0A6oHDgN2AO4BNwIUdql6SJElSwRg4cDcWLLiXI444kgMO+Giza374w6t57733mDPnLior+wKw664D2Xff/dlpp/7MmHElw4ePYPDg3Vm9+h1OPfV0jj56zNav/7//+1OeeebprAWsTO/BGga8FGNc0ejHRuBs4KkY4/djjK/EGC8GnkiPE0I4GBgFjI0xPhdj/CVwLjA5hNAz+ZckSZIkKR8dffQYhg0bzvTpl7N+/fom86tXr2bhwt9y/PEnbA1XDY0dO54ePUp58MH7APjyl4/j5JNPAWDTpk0sXPg7lix5jZEjD8raa2hLwHq5hblDgUcajT2SHt8yvyTGuLjRfAUwIsM/X5IkSVI3MHXqNFasWM7s2T9uMvfyy39j8+bNTe6z2qKsrIz99x/OCy88v834K6+8xGc+829Mm3Y+n/vc0RxySPb2cWjLJYK9QghPAXuSuh/rghjjn4BBwNJG65cBg9O/bmme9Jo/trHmbVRVVXTkcHUx9rPw2ePCZ48Lnz3uHuxzflm5spjS0rZtEN7W9R09LhPFxcUMGbInEyZ8g1tuuYnPfvaz7LPPfhQVFVFUVMR779UCsNNOO7ZYR79+/VixYtk284MHD2LOnDv5+99f4dprr2GnnXbijDMmtVpLe74PWg1YIYTewFCgmtSlfeuBScCjIYSPATsA7zc6bD2wZeuPJvMxxo0hhPoGa9qturqmo19CXURVVYX9LHD2uPDZ48Jnj7sH+5x/6urq2LSpLuP1paXFbVrfUHuPy8SW1/Gf/3kCv//977j88kuYNetO6uvrqa+vp7y8EoA1a2parKOmZg19+/bbZr5Pn0qGDq1k6NCP8Pbbq5gz5zZOPXUiJSUl262lpe+D4uIi+vcvb36utRcZY1wH7AgcHmN8PH3WahzwT+C/gXVA43upegJr079uMh9C6AEUNVgjSZIkSQCUlJQwdeo0Xn99CXPnzto6vu+++9OjRw+ef/7ZZo/buHEjL7304tZLCP/61z/z6qtxmzV77fVh1q9fz5o1a7JSe0bn92KMa2KM6xv8vg74G6lL/N4Adm10yEA+uCywpXloeumgJEmSJDF06F6MHTueefPmsGxZKjZUVlZy1FHHMH/+PN59d3WTY+bP/wnr1q1jzJgvA3DXXXO57bZbtlnz0kt/Y8cdd6Jfv35ZqbvVgBVC+HgIYU36csAtYyWkNqj4G7CI1PbrDR0OPJb+9SJgaAhhcKP5GqD56ClJkiSp2zvppHEMGTKUlSvf2jo2adLZ7LxzFRMnnsof/vA7VqxYzquvRq677mpmz/4xU6acx+677wHA8cefwJNP/h/z5/+EN998g1/84n7mz/8J48efTlFRUVZqzmSTi+eA14AfhxDOBGqB84CdgeuBXYA/hxAuBe4GTgAOAs5IH/8k8BRwbwhhUnr9lcC1McYOPamssk9ZRw6XJEmSROp9dSYPEe7s99+lpaVMnTqNCRPGbh3r06ecG2+cyc9+dg9z5tzG0qVvssMOffjoR0dw8823M2zYB4/uHTnyU3zve1cye/Zt3H77TD70oV341rfO5Zhjvpy1movq6+tbXRRC2A24Cvgs0Af4P2BKjPHF9PwX0vN7Aa8A58QYf9fg+AHALcBoUmeuZgMXpS81zMSewOJVq2qpq2u9XuUnb6gtfPa48NnjwmePuwf7nH9WrFjCgAF7ZLy+I5tcdBfb+zttsMnFEFIno7bKaJv2GONS4MTtzD8EPLSd+RXAVzL5syRJkiQpX2VvE3tJkiRJ6mYMWJIkSZKUEAOWJEmSJCXEgCVJkiRJCTFgSZIkSVJCDFiSJEmSlBADliRJkiQlxIAlSZIkSQnJ6EHDkiRJkgpX7byzqF+3ptV1Rb0rKT/5hqzUcNxxYygpKWHu3Hvo1avXNnOTJp3OoEGDOf/8iwF47721zJ8/j4ULf8uKFcspL69g2LDhnHDCyQwbNnybY3//+98wb94dvPnm6/TvvzPHHPNlTjjhZEpKSrLyOjyDJUmSJHVzmYSrtqxrr6VL32TmzJu3u2bNmneZOPEUfve73zB+/ETuuutnXHnltfTt25czz5zAQw89sHXtk0/+H5dddjFjxnyJuXPv4RvfmMRdd81l3rw5WXsNnsGSJEmS1CUMHLgbCxbcyxFHHMkBB3y02TU//OHVvPfee8yZcxeVlX0B2HXXgey77/7stFN/Zsy4kuHDRzB48O78/OcLOOywIzj22OMB2G23QSxZ8hoPPfQg48adlpXX4BksSZIkSV3C0UePYdiw4Uyffjnr169vMr969WoWLvwtxx9/wtZw1dDYsePp0aOUBx+8b+vvTzllwjZrioqKqKnJ3pk4A5YkSZKkLmPq1GmsWLGc2bN/3GTu5Zf/xubNm5vcZ7VFWVkZ++8/nBdeeB6AfffdnyFDhm6dX7u2lvvvX8BBBx2cneIxYEmSJEnqQgYP3p3x4ydyzz138sorL28zV1NTA0Dfvv1aPL5v376sXv1Ok/H333+fqVPPYf369ZxxxuRki27AgCVJkiSpSzn++BPZe+99uOKKS9m0adPW8b59U5cFrl1b2+KxtbU19Ou34zZjq1ev5pvf/G/+/vdXmDHjBgYM2DU7hWPAkiRJktTFlJSUMHXqNF5/fQlz587aOr7vvvvTo0cPnn/+2WaP27hxIy+99OI2lxAuX76Mb3zjVJYvX8pNN93Gvvvun9XaDViSJEmSupyhQ/di7NjxzJs3h2XLlgJQWVnJUUcdw/z583j33dVNjpk//yesW7eOMWO+DMA77/yLs876BvX1ddxyy2w+/OGPZL1uA5YkSZKkLumkk8YxZMhQVq58a+vYpElns/POVUyceCp/+MPvWLFiOa++GrnuuquZPfvHTJlyHrvvvgcAM2ZcyerVq7nkku/Ts2dPVq16m1Wr3uZf/1qVtZp9DpYkSZLUzRX1rszoIcJFvSs7oZoPlJaWMnXqNCZMGLt1rE+fcm68cSY/+9k9zJlzG0uXvskOO/Thox8dwc03386wYQcAsH79+zz22B+oq6vb5nhIXYL46KN/zErNRfX19Vn5wgnbE1i8alUtdXV5Ua/aoaqqgurqmlyXoSyyx4XPHhc+e9w92Of8s2LFEgYM2CPj9aWlxWzaVJfFivLf9v5Oi4uL6N+/HGAI8No2c1mvTJIkSZK6CQOWJEmSJCXEgCVJkiRJCTFgSZIkSVJCDFiSJEmSlBADliRJkiQlxIAlSZIkSQkxYEmSJElSQgxYkiRJkpQQA5YkSZIkJaQ01wVIkiRJyq3zF11GzYbaVtdVlJUzfdS0rNRw3HFjKCkpYe7ce+jVq9c2c5Mmnc6gQYM5//yLefPNN7j55ut4/vlnKSoqYsSIjzNp0rcYMGBAVupqK89gSZIkSd1cJuGqLevaa+nSN5k58+YW59etW8eUKZPYvLmO66+/lRkzbuLdd1dzzjlnsWHDhqzWlikDliRJkqQuYeDA3Viw4F5eeOG5Zuf/9KeneOutFXz3u5fz4Q9/hBD24aKLLuW11/7JSy+92MnVNs+AJUmSJKlLOProMQwbNpzp0y9n/fr1Teb3229/rrnmevr0Kd86VlycijQ1NWs6rc7tMWBJkiRJ6jKmTp3GihXLmT37x03mqqo+xMiRn9pm7M4776BXr14MHz6is0rcLgOWJEmSpC5j8ODdGT9+IvfccyevvPLydtfed9/PWLDgp5xxxmT69u3XSRVunwFLkiRJUpdy/PEnsvfe+3DFFZeyadOmZtfMnTuLGTOmc/LJp3Dsscd3coUtM2BJkiRJ6lJKSkqYOnUar7++hLlzZ20zV1dXx9VX/4DbbruFM86YzMSJZ+aoyuYZsCRJkiR1OUOH7sXYseOZN28Oy5Yt3Tp+7bVX8Ytf/JwLLvguJ544NocVNs+AJUmSJKlLOumkcQwZMpSVK98C4MknF3H//T/j618/lYMOOphVq97e+qO5XQdzwYAlSZIkdXMVZeWtL2rDuqSUlpYydeo0SkpKAPj1r38FwJw5t/GlLx21zY9HHvl9p9bWkqL6+vpc15CJPYHFq1bVUleXF/WqHaqqKqiursl1Gcoie1z47HHhs8fdg33OPytWLGHAgD0yXl9aWsymTXVZrCj/be/vtLi4iP79ywGGAK9tM5f1yiRJkiSpmzBgSZIkSVJCDFiSJEmSlBADliRJkiQlxIAlSZIkSQkxYEmSJElSQgxYkiRJkpQQA5YkSZIkJcSAJUmSJEkJKc11AZIkSZJy6x9TzmLzmjWtriuprGSva2/ISg3HHTeGkpIS5s69h169em0zN2nS6QwaNJjzz7+Y444bwzHHfIlx405r8jWOP/7LjB79ecaPnwjA888/yy233Mirr0bKyysYPfrzTJhwBj169MjKawDPYEmSJEndXibhqi3r2mvp0jeZOfPmRL7WihXL+fa3z2K//fZn7tx7uPDCS/j1r3/JrbfemMjXb4kBS5IkSVKXMHDgbixYcC8vvPBch7/W8uXLOOyww5k8eQq77TaIkSMP4jOf+SzPPPN0ApW2zIAlSZIkqUs4+ugxDBs2nOnTL2f9+vUd+loHHvhxLrro0q2/j/EVHn/8UT75yU91tMztMmBJkiRJ6jKmTp3GihXLmT37x4l9zaOO+nfGjz+JiooKxo0bn9jXbY6bXEiSJEnqMgYP3p3x4ycyc+bNHH74keyzz75N1sydO4u77prbZPz9999vMlZXV8cPf3gza9as4frrr+Gcc87mRz+6naKioqzUb8CSJEmS1KUcf/yJ/OEPv+eKKy5l1qw7m8z/x3/8J1/5yn82GT/77DOajBUXF7PvvvsDcOGFlzJx4jhefPF5Djjgo8kXjgFLkiRJUhdTUlLC1KnTGD/+JObOndVkvqKikkGDBjcZLy39IN4sXvxP3n57JSNHfnDP1V57fRiA6urqLFSd4j1YkiRJkrqcoUP3YuzY8cybN4dly5a2+fgnnnicSy65cJvNMl566UUA9txzSGJ1NmbAkiRJktQlnXTSOIYMGcrKlW+1+dijjvoCAFdccRlLlrzG008/xfTpl/OZz3yWoUP3SrrUrQxYkiRJUjdXUlmZ6LqklJaWMnXqNEpKStp8bP/+O3P99bfyzjvvcNppX+fyy7/Lpz99OBdeeGnrB3dAUX19fVb/gITsCSxetaqWurq8qFftUFVVQXV1Ta7LUBbZ48JnjwufPe4e7HP+WbFiCQMG7JHx+tLSYjZtqstiRflve3+nxcVF9O9fDjAEeG2buaxXJkmSJEndhAFLkiRJkhJiwJIkSZKkhBiwJEmSJCkhbXrQcAjhU8Ai4MgY4yPpsdHAVUAAXgXOizH+qsExHwJuAkYDG4A5wIUxxk1JvABJkiRJ6ioyPoMVQugDzANKGoztBzwA/A9wIPBz4P4Qwv4NDl0ADAAOA8YBpwDZ3RtRkiRJknKgLZcIXgu82WjsbOCpGOP3Y4yvxBgvBp5IjxNCOBgYBYyNMT4XY/wlcC4wOYTQs+PlS5IkSVLXkVHACiEcDXwBOKvR1KHAI43GHkmPb5lfEmNc3Gi+AhjRtlIlSZIkqWtr9R6sEMLOwO3AqcA7jaYHAUsbjS0DBrcyT3rNH9tSbPphXipgVVUVuS5BWWaPC589Lnz2uHuwz/ll5cpiSkvbtn9dW9d3N8XFxe36Pshkk4uZwIMxxodDCIMaze0AvN9obD3Qq6X5GOPGEEJ9gzUZW7Wqlrq6+rYepjzhU+MLnz0ufPa48Nnj7sE+55+6ujo2barLeH1paXGb1ndHdXV1LX4fFBcXtXjyZ7sBK4QwltTmFcNbWLIOaHwvVU9gbUvzIYQeQFGDNZIkSZJy6I4bn2Dd2o2truvdpwfjJh+SlRqOO24MJSUlzJ17D716bXsuZtKk0xk0aDDnn38xb775BjfffB3PP/8sRUVFjBjxcSZN+hYDBgzYur6+vp6HHnqABx64j8WL/0FJSQlDh36YL3/5OEaPPior9W/R2nnBcaQu81sRQqgFYnr8VyGEW4E3gF0bHTOQDy4LbGkeml46KEmSJCkHMglXbVnXXkuXvsnMmTe3/OevW8eUKZPYvLmO66+/lRkzbuLdd1dzzjlnsWHDBiAVri677GJuuOFajjjiSGbNupNbb53DIYeM4sorL+eKKy7L6mto7RLBk4DeDX4/AHgcOA34LfA9UtuvX95gzeHAY+lfLwKuDCEMjjG+0WC+Bni2Y6VLkiRJKiQDB+7GggX3csQRR3LAAR9tMv+nPz3FW2+tYM6cu+jTJ3WJ3kUXXcqxxx7DSy+9yIgRH+PBB+9n4cLfcvPNtzNs2AFbj91zzyGEsC9TpkziwAM/zlFHfSErr2G7ASvGuM1ZphDClvuplsYYV4YQbgT+HEK4FLgbOAE4CDgjve5J4Cng3hDCJGAX4Erg2hjjhuRehiRJkqR8d/TRY/jTn55i+vTLmT37Lnr23PZupP32259rrrl+a7iC1GYUADU1awBYsOCnHHLIqG3C1RYjRx7EyJEHsWDBvVkLWB3aOiTG+ALwFeA4UmekvgiMiTG+nJ6vT8+/RerM1xxgFpDd83KSJEmS8tLUqdNYsWI5s2f/uMlcVdWHGDnyU9uM3XnnHfTq1Yvhw0fw/vvv889//j+GDWtpCwk48MBP8MorL7NxY3Yud8xkF8GtYoxvktqgouHYQ8BD2zlmBamQJUmSJEnbNXjw7owfP5GZM2/m8MOPZJ999m1x7X33/YwFC37Kt751Ln379uPtt6upr6+nb9++LR7Tt29f6uvreffdd9l5550Tr9/N7yVJkiR1KccffyJ7770PV1xxKZs2bWp2zdy5s5gxYzonn3wKxx57PAAVFZUArF3b8obltbU1FBUVbTeEdYQBS5IkSVKXUlJSwtSp03j99SXMnTtrm7m6ujquvvoH3HbbLZxxxmQmTjxz61zPnj3ZZ5/9eP75lvfTe/bZvxLCvvTo0SMrtRuwJEmSJHU5Q4fuxdix45k3bw7Lln2w9961117FL37xcy644LuceOLYJsd99asnsGjRY7zwwnNN5p577q889dT/ceyxX81a3W26B0uSJEmSOstJJ43j0UcX8uqrfwfgyScXcf/9P+OUUyZw0EEHs2rV21vXlpdX0LNnT0aPPornn3+Wc845i/HjJ3LwwaO2Hnv77TM56qgv8PnPH5O1mg1YkiRJUjfXu0+PjB4i3LtPdi6ra0lpaSlTp05jwoTUmapf//pXAMyZcxtz5ty2zdqLL76Mz33uaADOOed8Djzw4/zv//50626Ee+31Ec49dyqf/exRWa25qL6+Pqt/QEL2BBavWlVLXV1e1Kt2qKqqoLq6JtdlKIvsceGzx4XPHncP9jn/rFixhAED9sh4fWlpMZs21WWxovy3vb/T4uIi+vcvBxgCvLbNXNYrkyRJkqRuwoAlSZIkSQkxYEmSJElSQgxYkiRJkpQQA5YkSZJUAPJk87q80JG/SwOWJEmSlOdKSkrZuHFDrssoGBs3bqCkpH1PtDJgSZIkSXmuvLwfq1dXs2HDes9kdUB9fT0bNqxn9epqysv7tetr+KBhSZIkKc/17t0HgHfffZvNmze1ur64uJi6Op+D1ZySklIqKnbc+nfaVgYsSZIkqQD07t0n41Dgw6Szx0sEJUmSJCkhBixJkiRJSogBS5IkSZISYsCSJEmSpIQYsCRJkiQpIQYsSZIkSUqIAUuSJEmSEmLAkiRJkqSEGLAkSZIkKSEGLEmSJElKiAFLkiRJkhJiwJIkSZKkhBiwJEmSJCkhBixJkiRJSogBS5IkSZISYsCSJEmSpIQYsCRJkiQpIQYsSZIkSUqIAUuSJEmSEmLAkiRJkqSEGLAkSZIkKSEGLEmSJElKiAFLkiRJkhJiwJIkSZKkhBiwJEmSJCkhBixJkiRJSogBS5IkSZISYsCSJEmSpIQYsCRJkiQpIQYsSZIkSUqIAUuSJEmSEmLAkiRJkqSEGLAkSZIkKSEGLEmSJElKiAFLkiRJkhJiwJIkSZKkhBiwJEmSJCkhBixJkiRJSogBS5IkSZISYsCSJEmSpIQYsCRJkiQpIQYsSZIkSUqIAUuSJEmSEmLAkiRJkqSElOa6AEmSJEnZc/6iy6jZUNvquoqycqaPmtYJFRU2z2BJkiRJBSyTcNWWddo+A5YkSZIkJcSAJUmSJEkJMWBJkiRJUkIMWJIkSZKUEAOWJEmSJCXEgCVJkiRJCcnoOVghhEHAD4HPkAplDwNTYozL0vMnAtOA3YHngMkxxqcbHP9h4CZgFPAOcEOM8eoEX4ckSZIk5VyrZ7BCCEXAQ8COwOHAYcCuwIPp+SOB2cAM4GPAC8BvQghV6fkyUoGsBvgkcB5wSQhhQtIvRpIkSZJyKZNLBHcBXgZOizE+F2N8DrgW+FgIYUfgXODuGOOPY4wvAxOBfwFbAtSxwADglBjjSzHG+cBVwDkJvxZJkiRJjVSUlSe6TtvX6iWCMcYVwNe2/D59ueBE4GngXeDfgEkN1teFEB4DDk0PHQo8E2Ns+GjoR0idxdolxvhWR1+EJEmSpOZNHzWtyVhVVQXV1TU5qKbwtWmTixDC/cAbwKeA04B+QB9gaaOly4DB6V8PamGeBmskSZIkKe9ltMlFA9OAHwAXAb8DRqbH32+0bj3QK/3rHYDqZuZpsCYj/ft72rLQVVVV5LoEZZk9Lnz2uPDZ4+7BPhc+e5wdbQpYMcbnAUIIXyN1Juuk9FTPRkt7AmvTv17XwjwN1mRk1apa6urq23KI8oinqgufPS589rjw2ePCc/6iy6jZUNvquoqy8mYvNVN+8nu5Y4qLi1o8+ZPJLoK7pAPVVjHG94B/AANJhaRdGx02kA8uC3yjhXloeumgJEmSOlEm4aot66TuLpN7sPYA7g4hfGLLQAihLxCAl4AnSG3dvmWuGPg08Fh6aBHwiRDCDg2+5uFAjDGu7Fj5kiRJktR1ZHKJ4DPA48DtIYTTgY3AdFL3Vc0ldSbrwRDCX4GFwBSgL3B7+vj7gO8D80MIFwEHkNra/cwEX4ckSZIk5VyrZ7BijHXAfwDPAr8AHgXWAIfFGGtjjA8DpwPfBv4C7AeMjjG+nT5+HXAUUElqa/fpwAUxxjsSfzWSJEmSlEMZbXKRDkvjtjM/B5iznfkIHNHW4iRJkiQpn7TpOViSJEmSpJYZsCRJkiQpIQYsSZIkSUqIAUuSJEmSEmLAkiRJ6sYqysoTXSd1dxntIihJkqTCNH3UtCZjVVUVVFfX5KAaKf95BkuSJEmSEmLAkiRJkqSEGLAkSZIkKSEGLEmSJElKiAFLkiRJkhJiwJIkSZKkhBiwJEmSJCkhBixJkiRJSogBS5IkSZISYsCSJEmSpIQYsCRJkiQpIQYsSZIkSUqIAUuSJEmSEmLAkiRJkqSEGLAkSZIkKSEGLEmSJElKiAFLkiRJkhJiwJIkSZKkhBiwJEmSJCkhBixJkiRJSogBS5IkSZISUprrAtQ9ffPGRaxZu6HVdZV9yrhu8qhOqEiSJEnqOAOWciKTcNWWdZKk5J2/6DJqNtS2uq6irJzpo6Z1QkWS1PV5iaAkSWpWJuGqLeskqTswYEmSJElSQgxYkiRJkpQQA5YkSZIkJcSAJUmSJEkJMWBJkiRJUkIMWJIkSZKUEAOWJEmSJCXEgKWcqOxTlug6SZIkqSsozXUB6p6umzyqyVhVVQXV1TU5qEaS1JyKsvKMHiJcUVbeCdVIUn4wYEmSpGZNHzWtyZgfhknS9nmJoCRJkiQlxIAlSZIkSQkxYEmSJElSQgxYkiRJkpQQA5YkSZIkJcSAJUmSJEkJMWBJkiRJUoHX6EgAABpxSURBVEIMWJIkSZKUEAOWJEmSJCXEgCVJkiRJCTFgSZIkSVJCDFiSJEmSlBADliRJkiQlxIAlSZIkSQkxYEmSJElSQgxYkiRJkpQQA5YkSZIkJcSAJUmSJEkJMWBJkiRJUkIMWJIkSZKUEAOWJEmSJCXEgCVJkiRJCTFgSZIkSVJCDFiSJEmSlBADliRJkiQlxIAlSZIkSQkxYEmSJElSQkozWRRC2AW4ChgN9Ab+CHw7xvhiev5EYBqwO/AcMDnG+HSD4z8M3ASMAt4BbogxXp3g65AkSZKknGv1DFYIoRi4D9gb+BJwCPAu8PsQQv8QwpHAbGAG8DHgBeA3IYSq9PFlwMNADfBJ4DzgkhDChORfjiRJkiTlTiZnsD4KHAzsF2N8GSCEcDLwL+ALwInA3THGH6fnJgJHABOAHwDHAgOAU2KMtcBLIYSPAOcAtyX7ciRJkiQpdzK5B+t14BggNhirA4qAHYF/Ax7ZMhFjrAMeAw5NDx0KPJMOV1s8AuydvvRQkiRJkgpCq2ewYoyrgIcaDZ8F9AKeAfoASxvNLwNGpn89qIV5gMHAW22oV5IkSZK6rIw2uWgohPBF4ArgWmBJevj9RsvWkwpgADsA1c3M02BNRvr3L2/LcuWhqqqKXJegLLPHhc8eFz573D3Y58Jnj7OjTQErhDCO1H1T9wDfIXWJIEDPRkt7AmvTv17XwjwN1mRk1apa6urq23KI8khVVQXV1TW5LkNZZI8Lnz0ufPa4e7DPhc8ed0xxcVGLJ38yfg5WCOFCYA5wK/D19L1W/yIVknZttHwgH1wW+EYL89D00kFJkiRJylsZBawQwneA7wHTYoyTY4z1AOmfnwAOa7C2GPg0qY0uABYBnwgh7NDgSx6eOjyu7PhLkCRJkqSuodVLBEMIw0lttz4buC2EMKDBdA2pe7EeDCH8FVgITAH6Aren19wHfB+YH0K4CDgAOBc4M6kXIUmSJEldQSZnsL4GlACnAssb/fhWjPFh4HTg28BfgP2A0THGtwFijOuAo4BK4GlgOnBBjPGORF+JJEmSJOVYJtu0XwBc0MqaOaTuz2ppPpJ6+LAkSZIkFayMN7mQJEmSJG2fAUuSJEmSEmLAkiRJkqSEGLAkSZIkKSEGLEmSJElKiAFLkiRJkhJiwJIkSZKkhBiwJEmSJCkhBixJkiRJSogBS5IkSZISYsCSJEmSpIQYsCRJkiQpIQYsSZIkSUqIAUuSJEmSEmLAkiRJkqSElOa6AElSfjp/0WXUbKhtdV1FWTnTR03rhIokSco9z2BJktolk3DVlnWSJBUCA5YkSZIkJcSAJUmSJEkJMWBJkiRJUkIMWJIkSZKUEAOWJEmSJCXEgCVJkiRJCTFgSZIkSVJCDFiSJEmSlBADliSpXSrKyhNdJ0lSISjNdQGSpPw0fdS0JmNVVRVUV9fkoBpJkroGz2BJkiRJUkIMWJIkSZKUEAOWJEmSJCXEgCVJkiRJCTFgSZIkSVJCDFiSJEmSlBADliRJkiQlxIAlSZIkSQkxYEmSJElSQgxYkiRJkpQQA5YkSZIkJcSAJUmSJEkJMWBJkiRJUkIMWJIkSZKUEAOWJEmSJCWkNNcFSCpMtfPOon7dmm3GappZV9S7kvKTb+icoiRJkrLMM1iSsqJxuOroOkmSpHxgwJIkSZKkhBiwJEmSJCkhBixJkiRJSogBS5IkSZISYsCSJEmSpIQYsCRJkiQpIQYsSZIkSUqIAUuSJEmSEmLAkpQVa+p6JbpOkiQpH5TmugBJheni1V/NeO3sLNYhSZLUmTyDJUmSJEkJMWBJkiRJUkIMWJIkSZKUEAOWJEmSJCXEgCUpKyr7lCW6TpIkKR+4i6CkrLhu8qgmY1VVFVRX1+SgGkmSpM7hGSxJkiRJSohnsJQTtfPOon7dmm3GmjuvUdS7kvKTb+icoiRJkqQO8gyWcqJxuOroOkmSJKkrMGBJkiRJUkIMWJIkSZKUEAOWJEmSJCXEgCVJkiRJCTFgSZIkSVJC2rxNewhhJlASYzytwdho4CogAK8C58UYf9Vg/kPATcBoYAMwB7gwxripY+VLknLFxy1IktRUxmewQghFIYTLgNMbje8HPAD8D3Ag8HPg/hDC/g2WLQAGAIcB44BTgEs7VLkkKad83IIkSU1lFLBCCEOBhcAZwOuNps8Gnooxfj/G+EqM8WLgifQ4IYSDgVHA2BjjczHGXwLnApNDCD0Teh2SJEmSlHOZnsE6GPgncACwuNHcocAjjcYeSY9vmV8SY1zcaL4CGJF5qSokRb0rE10nSZIkdQUZ3YMVY7wLuAsghNB4ehCwtNHYMmBwK/Ok1/wxw1rp378806Xq4qqmzMl1CcqRqqqKXJeghDR3v1VL7HthsZ/dg30ufPY4O9q8yUUzdgDebzS2HujV0nyMcWMIob7BmoysWlVLXV19e+tUF1dVVUF1dVvesinf2OPuy74XDr+Puwf7XPjscccUFxe1ePIniW3a1wGN76XqCaxtaT6E0AMoarBGkiRJkvJeEgHrDWDXRmMD+eCywJbmoemlg5IkSZKUt5IIWItIbb/e0OHAYw3mh4YQBjearwGeTeDPlyRJkqQuIYl7sG4E/hxCuBS4GzgBOIjUlu4ATwJPAfeGECYBuwBXAtfGGDck8OdLkiRJUpfQ4TNYMcYXgK8Ax5E6I/VFYEyM8eX0fH16/i3gcWAOMAu4rKN/tiQpd3zcgiRJTRXV1+fFrnx7AovdRbCwuZtN4bPHhc8eFz573D3Y58JnjzumwS6CQ4DXGs4lcYmgJKkb+uaNi1iztvUrvSv7lHHd5FGdUJEkSbmXxCYXkqRuKJNw1ZZ1kiQVAgOWJEmSJCXEgCVJkiRJCTFgSZIkSVJCDFiSJEmSlBADliRJkiQlxG3aJWXF+Ysuo2ZDbavrKsrKmT5qWidUJEmSlH2ewZKUFZmEq7askyRJygcGLEmSJElKiAFLktQulX3KEl0nSVIh8B4sSVK7XDd5VJOxqqoKqqtrclCNJEldg2ewJEmSJCkhBixJkiRJSogBS5IkSZISYsCSJEmSpIQYsCRJkiQpIe4iKCkrKsrKM3qIcEVZeSdUI0lS93XHjU+wbu3GVtf17tODcZMP6YSKCpsBS1JWTB81rcmYW3hLUtfjm+/Cl0l/27JO2+clgpIkSd2Yb76lZBmwJEmSJCkhBixJkiRJSoj3YEmSpGZ5b44ktZ1nsCRJUrO8N0eS2s6AJUmSJEkJMWBJkiRJUkIMWJIkSZKUEAOWJEmSVMB69+mR6Dptn7sISpIkdWO9+/TIeLdI5afmdvmsqqqguromB9UUPgOWJElSN+abbylZXiIoSZIkSQkxYEmSJElSQgxYkiSpWd4YL0lt5z1YkqR2uePGJzK+Mb65ezzU9XlvjiS1nWewJEntkkm4ass6SZIKgQFLkiRJkhLiJYKSJKlZtfPOon7dmm3Gmrs4sKh3JeUn39A5RUlSF+cZLEmS1KzG4aqj6ySpOzBgSZIkSVJCDFiSJEmSlBADliRJkiQlxIAlSZIkSQkxYEmSJElSQgxYkqR26d2nR6LrJEkqBD4HS5LULuMmH9JkrKqqgurq5p6UJKmr+uaNi1izdkOr6yr7lHHd5FGdUJGU3wxYkiRJ3Vgm4aot69T1/GPKWWxes+3z6v7ezLqSykr2utaHhneUlwhKkqRmranrleg6SbnROFx1dJ22zzNYkiSpWRev/mrGa2dnsQ5l16TFP6V88/utrqst6QUckf2CpDznGSxJkqRuLJNw1ZZ1UndnwJIkSZKkhBiwJEmSJCkh3oMlKSvuuPEJ1q3d2Oq63n16NLvdt6Tc894cSWo7A5Zy4vxFl1GzobbVdRVl5UwfNa0TKlLSMglXbVknqfN5b44ktZ2XCConMglXbVknSZIkdQUGLEmSJElKiAFLkiRJKmAllZWJrtP2eQ+WJKld/jHlLDavWbPN2N+bWVdSWcle197QOUVJkpq4cchXWbN2Q6vrKvuUcV0n1FPoPIMlSWqXxuGqo+sk5YZnNwpfJuGqLeu0fZ7BkiRJ6sZ2ORDq17W+rqh39muRCoFnsCRJkrqx+nWZnWXOdJ3U3RmwJEmSJCkhBixJWVGW4YNHM10nqfN5b44ktZ33YCknKsrKM3qIcEVZeSdUo2w4dPHdbVj9uazVIan9mtv9saqqgurqmhxUI6m9Lu/3UyqLW/9Ac01dL+CI7BdU4AxYyol9/nok69ZubHVd7z49YFQnFCRJklSgMglXbVmn7fMSQeVEJuGqLeskSZKkrsCAJUmSJEkJyatLBO+85Y+8+84HD2rosXkdp114VA4rkqTuq6SyMqOHCLsBgiSpO+m0gBVCKAG+B4wDKoCHgTNjjG+192tuLPGJd5KUK26AIElSU515ieAlwFjg68CngUHAgk788yVJktRIUe/MzjJnuk7q7jrlDFYIoQw4Gzgrxvjb9NjXgMUhhENijE90Rh2SJEnaVvnJno2WktRZZ7BGkLos8JEtAzHG14DXgEM7qQZJkiRJyqrOugdrUPrnpY3GlwGDO/KFq6oqOnK48oA9zk+L+/Vj4+rVra7r0a+fPS4w9rPw2ePuwT4Xjraci7TvHddZAWsHoC7G2PihRuuBXh35wp6+Lnz2OD8Nuea6JmMtXXJijwuHlxUVPnvcPdjnwlLUu5L6da3v+lrUu9K+Z6i4uIj+/cubneusgLUOKA4hlMYYNzUY7wms7aQa1IWUbXqPDaU7ZLROkiRJ7ed9dp2rswLWG+mfd23wa4CBNL1sUN3Av//rV214fs7ns1+QJEmSlIDOCljPkbr88zDgToAQwp7AnsBjnVSDuhCfnyNJkqRC1CkBK8a4PoTwI+CaEMLbwErgR8CjMcan2vt1e2xel1SJkiRJktRhnXUGC+AioAepM1g9gIeBM9vyBU464yDq6uqzUJokSZIkdVynBaz05hbfTv+QJEmSpILTWQ8aliRJkqSCZ8CSJEmSpIQYsCRJkiQpIQYsSZIkSUqIAUuSJEmSEmLAkiRJkqSEGLAkSZIkKSEGLEmSJElKiAFLkiRJkhJiwJIkSZKkhBiwJEmSJCkhpbkuIEMlAMXFRbmuQ1lmjwufPS589rjw2ePuwT4XPnvcfg3+7koazxXV19d3bjXtMwp4PNdFSJIkSVIDhwKLGg7kS8DqCYwElgObc1yLJEmSpO6tBNgVeBpY33AiXwKWJEmSJHV5bnIhSZIkSQkxYEmSJElSQgxYkiRJkpQQA5YkSZIkJcSAJUmSJEkJMWBJkiRJUkIMWJIkSZKUEAOWJEmSJCWkNNcFbE8IoQT4HjAOqAAeBs6MMb6Vy7rUfiGEQcAPgc+QCvgPA1NijMvS8ycC04DdgeeAyTHGp3NUrtophHAa8B1gMPAScG6McWF6bjRwFRCAV4HzYoy/ylWtap8QQiWpPn4R6AX8ktT38sr0vH3OYyGEmUBJjPG0BmOTgEmkvq+XANfGGG9vMP9h4CZgFPAOcEOM8epOLVwZa6HHTwOfaLR01pY1IYQPkerxaGADMAe4MMa4qXOqVlu00OP/JPU+ayip7+OrY4xzGszb4wR09TNYlwBjga8DnwYGAQtyWZDaL4RQBDwE7AgcDhwG7Ao8mJ4/EpgNzAA+BrwA/CaEUJWTgtUuIYSxwM3AdOAA4FHggRDCniGE/YAHgP8BDgR+DtwfQtg/V/Wq3f4H+DxwCnAoUA78IYTQ0z7nrxBCUQjhMuD0RuNnkPqe/h4wHLgW+FEI4eT0fBmpD8xqgE8C5wGXhBAmdGL5ysB2elwE7AucSOr/zVt+TGmwbAEwgNT/v8eR+v6/NPtVqy220+NDgbtIBagDgOuB20IIX2iwzB4noKi+vj7XNTQr/Y/128BZMcY70mN7AouBf4sxPpG76tQeIYQBwHXA+THG19JjXwLuB3YC7gGWxxjHpeeKSX3yPSvG+INc1Ky2Sf8PejHwkxjjtPRYMfAXUmczDgNCjPHfGxzzB+DVGOPpTb+iuqIQwgjgr8BnY4y/S4+VA28A3wQOwT7nnRDCUGAWMAx4D/htgzMXzwEPxxjPa7B+FjAkxnhECOG/gNuAATHG2vT8d4ETYoyhk1+KWtBKj/cC/h8wNMa4uJljDwaeaDif/kDtRqAqxri+c16FtqeVHl8DHB5j/HiD9X8CnooxnmWPk9OVz2CNIHVZ4CNbBtJvyl8j9Wmp8kyMcUWM8WsNwtUgYCLwNPAu8G9s2+864DHsdz4JwB7AvVsGYox1McYRMcb5pHr5SKNjHsEe55uPpH9etGUg/ab6VVIh2j7np4OBf5L6ZLvxG+yzgFsbjdWRuiIBUr19Zku4SnsE2DuEsEvypaqdttfjYcA6UpeNNedQYEmj8PUIqfdqI5ItUx2wvR5XA/uHEA5Pn+X6NKm+P5Oet8cJ6cr3YA1K/7y00fgyUtd/K4+FEO4HvkTqOv1/B/oBfWi+3yM7tTh1xN7pn/uFEBaS+of7FVJnLZ8g9X3t93T+W5b+eRCpT7y33DM7CFiJfc5LMca7SF0+RAih8dyjDX8fQtgd+C9Sn2xDyz2HVN+9d7oL2F6PSf17vRq4K4RwGLCK1P0316U/8Gytx3/MUtlqg1Z6fDOpD7MXApuBEuCaGONP0vP2OCFd+QzWDkBdjHFjo/H1pG6oVn6bBhxE6hPw35H6dATg/Ubr7Hd+qUz/PBe4HTgKeBFYGELYl9T3tT3Of0+TCs63hhB2DSH0Bq4AqoAy7HNBS98X+xCwgtR9WdByz8G+54v9Sd1L+Wvgc6TejF8KfDc936TH6fdo9djjfPEhYBdSm1B9gtSZ6TNDCKem5+1xQrryGax1QHEIobTRziU9gbU5qkkJiTE+DxBC+Bqp+zZOSk/1bLTUfueXLR+IfD99SSAhhDNJXXZwBqnva3uc52KMG0IIXwHuJPXp5gZSn5j+Mv1r+1yg0vd3/IrUG7HDYozvpqda6jnY93zxdaA8xrg6/fsXQgh9gQtDCJfQTI9DCD2AIuxxvrgN+GuD3T2fTX9gclUIYQ72ODFd+QzWG+mfd200PpCmpy+VB0IIu6QD1VYxxveAf5Dq61rsd77b0qsXtgzEGOuBl4EhpL6v7XEBiDG+EmP8BLAzsHOM8VRSl5D8A/tckEIIB5K6Ab4OOCTG+M8G0y31HOx7XogxbmoQrrZ4gdQVJn2xx4XgU3xwv9UWfwT6k7pVwx4npCsHrOdIbfd62JaB9C6Ce5La+ED5Zw/g7hDC1mdspD8dC6SelfQE2/a7mNT2/PY7f/yFVFDeet9cemfB/Ui98V5Egx6nHY49zishhIoQwqMhhGExxlUxxpr0v88fBX6DfS44IYR9SF3OvQQYFWN8o9GSRcAnQgg7NBg7HIhbno2mri2E8FQI4bpGw58AlqWD1yJgaAih4b2Uh5N6r/ZsJ5WpjnmT1GMWGhoGrIoxvoM9TkyXvUQwxrg+hPAj4JoQwtukbpz+EfBojPGp3FandnoGeBy4PYRwOqnLyaaT2tVmLqk34A+GEP5K6gbMKaQ+Nbu9+S+nribG+F4I4YfA90MIb5H69PO/gb2AY0ndn/PnEMKlwN3ACaTuxTsjRyWrHdKBqgS4LoRwFqn7NmYDv48xLgwhVGOfC81PSN2bcTLQI/3YDYBNMca3gfuA7wPzQwgXkdrB7FzgzFwUq3b5X+CyEMJfgP8jtQHVecDZ6fkngaeAe9MPnd4FuJLUA6c3dH65aofrgR+GEF4ida/dwcAFwGXpeXuckK58BgvgIlLX9d8J/IHUJ2fH5bQitVt6F6L/IPUpyC9IPYB2Danr+GtjjA+Teijet0mdCdkPGJ3+n7fyxzTgalLPPHuB1D/go2PKC8BXSH0fPwt8ERgTY3w5V8Wq3b4G1JL6H/IDpL6f/wPAPheWEMLepM5KDwQisLzBj6cAYozrSG1qU0lqE5TpwAVbnmOpvHA1qTfbFwF/IxWuvhVjvB22Xu79FVI7Qj5OaofBWXzw5lxdXPz/7d1PqKVzGAfwL6kpWfhTslC3IT3RJEVKkaJEUsrsKLPQsLAwRKEYyYq5iAillI2wUFaymA0LpTRFnoiZnQWZBWOmMV2L9526XXMv8s7cc858PpvnnOc95/S8ncXpe37nPb/u1zJ86fFghj+gejLDe748HvceT2RmNxoGAACYN7O+ggUAADA3BCwAAICJCFgAAAATEbAAAAAmImABAABMRMACAACYyMxuNAzA6a2qdid5+l8+/ECS3Rn2bdnV3S+dpLEAYEMCFgCzau8JejuSLCV5OcnBVf2DGTY1fibj5rcAsBlsNAzA3KiqvUluTLK1u/dv7jQA8HeuwQIAAJiInwgCsBCqakfWXINVVfuTfJ9kV5Lnk1yf5HCSD8feeUmWk9ya5EiST5I81N0/r3ntm5I8nuTaDJ+d+5Ls6e4PTvJpATBnrGABsOi2Jvksw2fe60l+SrIzyTtjfynJmxmC2N1J3lr95Kq6L8mnSa5M8l6SN5JcmOT9qnri1JwCAPNCwAJg0V2S5O3uvqW7H82winUoyfYkXyS5buzfkCFk3VlVZydJVV2c5NUk3ya5ort3dvfDSbYl+TzJs1W17ZSfEQAzS8AC4HTw4vEb3X0wyTfj3eXuXhn7x5J8OfaXxnpPki1JnuruX1a9xh8Z/kL+zCT3ntzRAZgnrsECYNEd7e4Da3q/j/XHNf3DY90y1qvHevMJVqrOGetV/39EABaFgAXAoju0wbEj//Dcc8f6wAaPOf+/jQPAIhOwAGB9v4310u7+YVMnAWAuuAYLANa3b6zXrD1QVZdV1QtVdccpngmAGSZgAcD63k1yLMlzVXXR8WZVnZXklSSPJLlgk2YDYAb5iSAArKO7v6uqx5LsSfJ1VX2U5NcktyW5PMnHGUIYACSxggUAG+ru5SS3J/kqyV1J7k9yNMPq1fbu/nMTxwNgxpyxsrKy2TMAAAAsBCtYAAAAExGwAAAAJiJgAQAATETAAgAAmIiABQAAMBEBCwAAYCICFgAAwEQELAAAgIkIWAAAABP5Cy8qvrVphXuXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "tickparams = {'fontsize':15, 'xticks':catalysis_data.Time.values}\n",
    "markerparams = {'marker':'s', 'markersize':10}\n",
    "ax=catalysis_data.plot(style='s', x=0, figsize=(12, 7), xlim=(0, 200),\n",
    "                       **tickparams, \n",
    "                       **markerparams)\n",
    "ax.set_xlabel('Time', fontsize=20)\n",
    "ax.legend(loc='best', fontsize=15)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The observed data is the flattened vector of concentrations of each of the reactants and products at $t=20, 60, 90, 120, 150, 180$. We collect the observed data and flatten the data matrix into a vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y = catalysis_data.values\n",
    "data = Y[1:, 1:].flatten()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also require the model (i.e. solver) for the dynamical system governing the catalytic conversion. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def A(x):\n",
    "    \"\"\"\n",
    "    Assemble the reaction coefficient matrix. \n",
    "    \"\"\"\n",
    "    k = tt.exp(x) / 180.\n",
    "    res = tt.zeros((6, 6))\n",
    "    res = tt.set_subtensor(res[0, 0], -k[0])\n",
    "    res = tt.set_subtensor(res[1, 0], k[0])\n",
    "    res = tt.set_subtensor(res[1, 1], -(k[1] + k[3] + k[4]))\n",
    "    res = tt.set_subtensor(res[2, 1],  k[1])\n",
    "    res = tt.set_subtensor(res[2, 2],  -k[2])\n",
    "    res = tt.set_subtensor(res[3, 2],  k[2])\n",
    "    res = tt.set_subtensor(res[4, 1],  k[4])\n",
    "    res = tt.set_subtensor(res[5, 1],  k[3])\n",
    "    return res\n",
    "\n",
    "def g(x, z):\n",
    "    \"\"\"\n",
    "    The right hand side of the dynamical system.\n",
    "    \"\"\"\n",
    "    return tt.dot(A(x), z)\n",
    "\n",
    "def euler(x, zn, dt):\n",
    "    return zn + dt*g(x, zn)\n",
    "\n",
    "def Z(z0, x, T):\n",
    "    \"\"\"\n",
    "    Return the solution of the dynamical system for a given time discretization. \n",
    "    \n",
    "    We solve the ODE by applying a simple first order Euler discretization. \n",
    "    \n",
    "    INPUTS:\n",
    "        1. z0 -> Initial conditions.\n",
    "        2. x  -> Unscaled rate coeff. vector.\n",
    "        3. T -> Time discretization. \n",
    "    \"\"\"\n",
    "    N  = len(T)  # number of time steps \n",
    "    t0 = T[0]\n",
    "    dt = T[1]-T[0]\n",
    "    fn = lambda z : euler(x, z, dt)  # update scheme of discretized ODE\n",
    "    values, _ = th.scan(fn=fn, outputs_info=z0, non_sequences=[], n_steps=N-1)\n",
    "    return values\n",
    "\n",
    "def CatalyticSolver(x):\n",
    "    \"\"\"\n",
    "    Obtain a prediction of the reactant/product concentrations as a function \n",
    "    of the reaction rate coefficients. \n",
    "    \n",
    "    INPUTS:\n",
    "        x (theano.tensor) : 6-dimensional array of reaction rate coefficients. \n",
    "    \n",
    "    RETURNS:\n",
    "        The observation of the concentrations for 5 reactants/products at \n",
    "        6 observation times (30 dimensional vector). \n",
    "        The order of the output matches the order of output in the observed dataset.\n",
    "    \"\"\"\n",
    "    z0 = np.array([500.,   0.,   0.,   0.,   0.,   0.])\n",
    "    z0_s = shared(z0)\n",
    "    T = np.arange(0., 180.01, 0.5)\n",
    "    Tobs = np.array([30.*i for i in range(1, 7)])\n",
    "    tobsidx = [np.where(T[1:] == tobs)[0][0] for tobs in Tobs]\n",
    "    zidx = [0,1,3,4,5]\n",
    "    Zp_all = Z(z0_s, x, T)  # entire simulation result\n",
    "    Zp_pred = Zp_all[tobsidx][:,zidx]\n",
    "    Zp_pred_flat = Zp_pred.flatten()  # flatten the result to match observed data shape\n",
    "    return Zp_pred_flat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that we are working with the transformed variable `x` which is converted into the reaction rate vector `k` inside the solver. \n",
    "\n",
    "Suppose we have the following Gaussian prior for the reaction rates:\n",
    "$$\n",
    "p(\\mathbf{x}) = \\mathcal{N}(\\mathbf{x}|\\mathbf{0}, \\gamma^2 \\mathbf{I}),\n",
    "$$\n",
    "and the measurement process is defined by the following likelihood model:\n",
    "$$\n",
    "p(y|\\mathbf{x}, \\sigma^2, \\gamma) = \\mathcal{N}(y|f(\\mathbf{x}), \\sigma^2),\n",
    "$$\n",
    "where, $f(\\cdot)$ denotes the dynamical system model. \n",
    "The likelihood specification encodes the assumption that observations (i.e. concentrations of all substances at all observation times) are independent of each other and individually follow a Gaussian distribution with fixed variance. \n",
    "\n",
    "The likelihood mean $f(\\cdot)$ has already been implemented for you as the `CatalyticSolver` function. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part A -  Constant $\\gamma$ and $\\sigma$\n",
    "\n",
    "Fix the prior precision parameter $\\gamma$ and the likelihood noise $\\sigma$ as constants and use `PyMC3` to estimate the posterior distribution over the transformed reaction rates vector $x$ and use the posterior distribution to show the predictive distribution of reactact/product concentrations as a function of time.  \n",
    "\n",
    "Here is the workflow you should follow:\n",
    "\n",
    "1. Express the model using suitable distributions for the priors and the likelihood with a `pymc3.model` context.\n",
    "\n",
    "2. Simulate a Markov Chain to generate samples of $x$ (you just need to call `pymc3.sample` with appropriate arguments). \n",
    "\n",
    "3. Perform model diagnostics, i.e., tune the hyperparameters, $\\gamma$ and $\\sigma$, burn early samples of the chain to remove transient phase samples, and thin the MC to remove autocorrelation.\n",
    "\n",
    "4. Plot the posterior predictive distribution of the reactant/product concentrations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution:**\n",
    "<br></br><br></br><br></br><br></br><br></br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part B -  Exponential prior over $\\sigma^2$ and $\\gamma^2$\n",
    "\n",
    "We will now introduce prior specifications for $\\sigma^2$ and $\\gamma^2$. Let $p(\\sigma^2) = \\mathrm{Exp}(\\sigma^2 | \\alpha_1)$ and $p(\\gamma^2) = \\mathrm{Exp}(\\gamma^2 | \\alpha_2)$, where, $\\mathrm{Exp}$ denotes the exponential distribution. This formulation introduces 2 additional parameters to tune - $\\alpha_1$ and $\\alpha_2$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution:**\n",
    "<br></br><br></br><br></br><br></br><br></br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part C - A model with a different noise variance for each observed species with Jeffrey's prior."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now consider the following likelihood:\n",
    "$$\n",
    "p(y_i|\\mathbf{x}) = \\mathcal{N}(y_i|f(\\mathbf{x}), \\sigma_{i}^{2}),\n",
    "$$\n",
    "where $i$ is an index for observed chemical species.  In other words, we construct the likelihood model such that there is a different noise parameter, $\\sigma_i$ associated with the measurements obtained for each different species. \n",
    "Each $\\sigma_{i}^{2}$ is specified with a Jeffreys' prior, i.e., $p(\\sigma_{i}^{2}) \\propto \\frac{1}{\\sigma_{i}^{2}}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution:**\n",
    "<br></br><br></br><br></br><br></br><br></br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part D - Model Selection using SMC\n",
    "\n",
    "Use Sequential Monte Carlo (SMC) to determine the model evidence of the 3 different models defined in parts A to C. Which is the best model that you find? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution:**\n",
    "<br></br><br></br><br></br><br></br><br></br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 2  - Bayesian Linear regression\n",
    "\n",
    "Recall that the standard Bayesian Linear regression model admits closed form expressions for the posterior distribution over the weights and the posterior predictive distribution over the observations. Suppose you select a vector of suitable basis functions $\\phi(x) = (\\phi_1(x), \\phi_2(x) \\dots, \\phi_M(x))^T$, the predictive model is the GLM:\n",
    "$$\n",
    "f(x) = \\mathbf{w}^T \\phi(x).\n",
    "$$\n",
    "\n",
    "If weights are equipped with a Gaussian prior $p(\\mathbf{w}) \\sim \\mathcal{N}(\\mathbf{w}|0, \\Sigma_p)$ and the observations are assumed to follow a Gaussian distribution - $\\mathbf{y} \\sim \\mathcal{N}(\\mathbf{y} | \\Phi \\mathbf{w}, \\sigma^2 I)$, the posterior distribution over the weights is given by - $p(\\mathbf{w} | \\mathbf{X}, \\mathbf{y}) = \\mathcal{N}(\\mathbf{w} | \\mathbf{m}, \\mathbf{S})$, where, $\\mathbf{S} = (\\sigma^{-2} \\Phi^T \\Phi + \\Sigma_p )^{-1}$ and $\\mathbf{m} = \\sigma^{-2} \\mathbf{S} \\Phi^T \\mathbf{y}$ and the posterior predictive distribution over the observations are given by $p(y^*|x^*, \\mathbf{X}, \\mathbf{y}) =  \\mathcal{N} ( y^* | \\mathbf{m}^T \\phi(x^*), \\phi(x^*)^T \\mathbf{S} \\phi(x^*) + \\sigma^2)$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part A - Compare MCMC to analytical solution\n",
    "\n",
    "Setup a GLM for the motorcycle data (loaded below), with an fixed precision prior on the weights and a constant likelihood noise, and use `PyMC3` to determine the posterior over the weights and the posterior predictive distribution at new test inputs. Compare your MCMC solution for the posterior with the analytical solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x139d50748>]"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD7CAYAAACIYvgKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAYPklEQVR4nO3dfZRcZ33Y8e/sroyRtZGl9fhYrBCLLfwoyZFjYoxx/FI3cUk4bks4vDg2L6anCaENNGkLtKfAgTgxpycQkzYlCSeQQ3ugNseUpAkKOY7BBmPJkNNjRzKOHizZG0XrdSuvFMWysdHuTP+YGWm0zO7Oy71zX+b7+Uea0Z2599G98/zu8/ye57mVer2OJGm0jWV9AJKk7BkMJEkGA0mSwUCShMFAkgRMZH0AXXgRcDkwDyxlfCySVBTjwBbgr4AX1tq4CMHgcuD+rA9CkgrqGuBba21UhGAwD3Ds2LNs2nQOCwsnsj6e1ExNbSht+SxbMVm2Ypqa2nCqzqRZh66lCMFgCaBWa0yOa/1ZVmUun2UrJstWTG1l66p73QSyJMlgIEkyGEiSMBhIkjAYSJIwGEiJOzB3nF17ZjkwdzzX3ym1K8LQUqkwDswd5+N3PMTiUo2J8THef9Mr2T69MXffKS1ny0BKUDx0jMWlGvU6LC3ViIeO5fI7peUMBlKCwrZNTIyPMVaB8fExwrZNufxOaTm7iaQEbZ/eyPtveiXx0DHCtk2JdOek8Z3ScgYDKWHbpzcmXmEv/84Dc8cNDkqUwUAqGBPKSoM5A6lgTCgrDQYDqWBMKCsNdhNJBWNCWWkwGEgFlEaSWqPNbiJJksFAkmQwkCRhMJAkYTCQJGEwkCRhMJAkYTCQJGEwkIbGR1cqz5yBLA2BK40q72wZSEPgSqPKu1RaBiGEe4HzgZPNt34ZuAj4ELAO+J0Y46fS2LeUR62VRpeWaq40qlxKPBiEECrAxcDLYoyLzfemgTuBy4AXgN0hhHtjjI8mvX8pj1xpVHmXRssgNP+8O4QwBfwh8Azw9RjjUYAQwpeANwG3prB/KZdcaVR5lkYw2AR8DXgvjS6h+4AvAvNt28wDr+7lS6emNgBQrU4mcYy5VebyWbZismzF1Kozu5V4MIgx7gH2tF6HED4L3A78ZttmFaDWy/cuLJxgamoDR448k8hx5lG1Olna8hW1bN08eL7fsmXxUPte91nU89aNspetVWd2K42cwdXAi2KMX2u+VQFmgS1tm10APJn0vqUkpTkcNIuhpg5v1WrSGFp6LvDxEMLZIYRJ4BbgbcDPhBCqIYT1wBuBv0hh31Ji0hwOutp3pzU5zeGtWk0a3URfCSFcATwEjAOfijE+EEL4IHAvcBbwmRjjd5Let5SkNIeDrvTdad69O7xVq6nU6/Wsj2EtM8AT5gyKrahlG3bOYNeeWb78zcep12GsAm+49kJuuHJmgBKsvc/VFPW8daPsZWvLGbycRlf9qlyOQlpFmsNBO3132nfvDm/VSgwGUo44OU1ZMRhIOePdu7LgQnWSJIOBJMlgIA1s/+xRH1qjwjNnIA3gwNxxPnHnQ5xcdFavis2WgTSAeOgYi4vZzOr1MZpKki0DaQBh2yYmJsZYXBzurF7XGVLSDAbSALZPb+S2d1/Fg3vnhjovoNM6QwYDDcJgIA1ox8xmps5ZN9R9us6QkmYwkAooLzOVs3gmg9JhMJAKKuuZyp3yFoDBoaAMBpL6sjxvsXvfPA888pRJ7YJyaKmkvrTyFmMVGB9vVCU+PKe4bBlI6svyvAXAA488ZVK7oAwGkvq2PG+Rh6S2+mMwkJSYrJPa6p85A0mSwUCSZDCQJGEwkCRhMJAkYTCQJGEwkCRhMJBywyeXKUtOOpNyoKxPLnOJ6+IwGEg5UMYnl5U1wJWV3URSB8Puslm+AmgZFnnrFOCUX7YMpGWyuKPNy5PLkuSjOYvFYCAtk1WXTdkWees1wK2VX8gq/5Bl3mOY+zYYSMt4R5ucbgPcWq2xrPIPWeY9hr1vcwbSMq072jdce2HfP0CHifZmrfxCVvmHLPMew963LQOpg0G6bBxF07u1WmNZtdaybCUOe98GAylhZRwmmra18gtZJdizTOwPe98GAylhRc05ZD1BbK3WWFYJ9iwT+8Pct8FAakqqMiziMFG7tmQwkEi+MizaMFG7tjTU0UQhhJtDCI+GEB4LIfzKMPctrSavs2WHNSqpjDOg1ZuhtQxCCNPAbcBlwAvA7hDCvTHGR4d1DNJK8tjPP8yumyJ2bSlZw+wmuh74eozxKEAI4UvAm4Bb09zpSv3AWSfLetHvsRapjFnLY2U47K6bJLu2urn28nJ95uU4sjbMYPASYL7t9Tzw6jR32LqzOrlYY2yswtteezHXXTqd6B1XLxdSPxddv8dqQrB3eevnz6K1kkTF2M1s4t375vnWvnmWavVMr09/J6cNMxiMAfW21xWg1u2Hp6Y2AFCtTna9w/v2zrO42NhFrVbnC3d/j8kNZ7N775On7rgWF2t89duHuPlnd7BjZnNX37t/9ij7Dj7Ns98/yZ984yC1ep11E2Pc9u6r2DGz+dS/77zoPAD2HXyayfVn8Yf/ex+LizUmmtu2/q213X17v3fGZ3ZedB6HF55jqe3u8PDCc1x56dZVj2ty/VlnlLH9c+3H1m15k9LLuSuaNMpWrU7ysXPXD+187Z89yifufOiMa7Ra7b1s9+2dX/Gabe3j5Mnaqcpgres6Tb38vgbR/tt85rkfDOV8turMbg0zGBwGrml7fQHwZLcfXlg4wdTUBo4ceabrHW6dWk9lrEK91rjslmp1/uDLe1mqnY5JdeDh7x3hkccXuroraG9ttDu5WOPBvXMc+/vnTt1pjI1VqDT3W6k0jqNOIwDtuv8gDzzy1Bnb1Wp1Km2fmRgf46brX8H4+Bg07w63Tq3v+H+w0nFVmgnBrVPr2fPw4czugqrVyZ7OXZGkWbapc9Zx3SVbAFL//3tw7xwnF0/fJD24d44dM5t73u/WqfUrXrOn9tG2/WrXdZqq1clVjzUpy3+bFWBiIt3fX7U6earO7NYwg8E9wEdDCFXgWeCNwLvS3OH26Y287bUX84W7v0et3qiQa81AUAGqm87myLHnqdN9n2yrH3e5sUqFsG3TGf28taX66Yu+XmdsrEK9Xm9cfNBxu0rb35eWajz7/ZNd9WV3Oq4K8GMzm3j91ReyfXoju/bMOnxQK0qqW2q1/Ev7PipjFa7ZuYWf2rkl00ldaeeKlv82e6lvhmlowSDGOBdC+CBwL3AW8JkY43fS3u91l06ztbqBeOgY57x4HXfc89ipi/3nrnjZGa+7ufhbF/Ni293N2FiFt7724lMntv1ib93xjzfv8p/9/slT+3ngkafO2G5pqQ4VGK80gkZlrMLC8ecBuOHKmZ6Oq1JpHEcrELRvs1Z5y5hQK2OZkpZkxbhS/iWPifq0c0Wdfpt5GbHWrlKv19feKlszwBP9dBN1srxS6Dep2wourcp9pZFKwIrf377d4SMnGi2YWp2x8QqXXDjFvscXekqwrXVcncrf6TvS6ErKspso7SShXWDFNMyydfPbTNKybqKXA7NrfWbkZiAvvwvo566g1zVUVtq2fbt46Bj1eqOLqF6r84PFJZZq9Z66dLopy1rblHEmahnLpGLJ20i1TnyeQU6EbZuYmDg9A/SycH4mM0LLOBO1jGWSkjZyLYO82j69kdvefRUP7p071YRs5TpGZcnetJSxTFLSDAY5smNmM1PnrDv1ehSX7E1LGcskJcluIkmSwWCU+Fxe9aPI102Rj33Y7CYaEa7Bon7snz1a2OvGa743tgxGRF7X61e+7Tv4dGGvm2Fc82VqedgyGBF5XK8/a0WdlTzM49550XmFvW7SvubL1vIwGIwIh1eeaZAfcpZBZNgV0I6ZzYW9btK+5ss2mdFgMEIcXnlavz/kTpXxMJfmzqICKvJ1k+axl621bTDQSOr3h9ypMh7mOvxlq4CKrGytbYOBRlK/P+SsK+OyVUBFV+RW03IGAwHDX1UxD9b6IXfKDeShMh5WBXRg7jj37Z1n69T60l8LMhiMtPYAcMc9jw31SUx5t1qitkx3gytplb/VAhrla2FUGAxGVHtl13okZ0ten8Q0TGUbKdKr9vIzguUfRU46G1HtP/Z6rX7qOcyQ3ycxDdOoL3s96uUfRbYMRtTyRGjrkZyjlDNYTR5yA1lqlf/wwnPmDEaEwWBEjUplN8gEsVHIDaxm+/RGrrx0ayEee1nU2eR5YjAYYWWv7Mq2XIA68zwnw5yBSsvF+UaD5zkZBgOVlknQ0eB5TobdRCqtUcmLjDrPczIMBiq1sudF1OB5HpzdRJIkg4EkyWAgScJgIEnCYKAelOnh32XjudGgHE2krjjLM788N0qCLQN1xVme+eW5URIMBiOu2+4FZ3nml+dGSbCbaIT10r3gLM/88twoCQaDEba8e2H3vvlVKxRneeaX50aDMhiMsPYH3FTGKnxr3zxLtbpJSGkEmTMYYa3uhTdceyHX7NzCUq1uElIaUbYMRlyre+HA3HEeeOSpU4/BHJUkpE/IkhoMBgJGMwnp+HzptMSDQQjhFuA/A/+3+dauGOMHQwjbgM8D5wMReGuM8UTS+1f/Ri0J2Wl8/iiVX2qXRsvgVcC/izHesez93wN+L8Z4Zwjhw8CHgf+Qwv6lrrQn0Eepa0zqJI1gcDnwihDCfwL+GngvcAK4Fvj55jafA76BwUAZGsWuMWklaYwmmgd+A7gE+DvgvwHnAf8QY1xs22ZrCvuWerJ9eiM3XDljINDIq9Tr9b4+GEJ4M/DJZW/vjzFe37bNJuAgsBN4MMb40ub7E8CJGOPZXexqBniir4OUJL0cmF1ro767iWKMdwF3tb8XQtgYQvi3McZWkKgAi8D/AzaGEMZjjEvAFuDJXva3sHCCqakNHDnyTL+HnHvV6mRpy5d22bIcIup5K6ayl61VZ3Yr6W6iE8AHQghXNF+/B/jjGONJ4H7gxub77wC+mvC+NaJaQ0S//M3H+fgdD7mmv9SHRINB867/LcDvhxD+BrgM+EDzn/818K4QwqPANcCHkty3RpdLOEuDS3w0UYzxfuAnO7z/t8B1Se9PcoioNDhnIKvwHCIqDc5goI6KtmbPqM2elpJmMNAPcc0eafS4hLV+iAlZafQYDPRDfKauNHrsJtIPMSErjR6DgToyISuNFruJVCoH5o6za8+ss5ClHtkyUGl0MwqqaENmpWExGKg01npymUNmpZXZTaTSWGsUlENmpZXZMlBprDUKyjWMpJUZDFQqq42CcsistDKDgVZUxmSrQ2alzgwG6qjoydYyBjIpTQYDdbTWyJw8K3ogk7LgaCJ1VOT1iRw1JPXOloE6KnKy1VFDUu8MBlpRUZOtRQ5kUlYMBiqlogYyKSvmDCRJBgNJksFAkoTBQJKEwUCShMFAkoTBQJKEwUAp8DnEUvE46UyJcpE4qZhsGShRLhInFZPBQIkq8mqn0iizm0iJcpE4qZgMBkqci8RJxWM3kSTJYCBJMhhIkjAYqOCc4CYlwwSyCssJblJybBmosJzgJiVn4JZBCOE3gKUY40ebr88FvgBcCBwB3hJjfCqEcBbwWeBVwPeBm2OM+wfdv4bvwNzxXMwjaE1wW1qqOcFNGlDfwSCEsBG4HbgJ+K22f/pN4P4Y4w0hhLcD/wW4Efg3wLMxxh8NIVwLfA54Tb/7Vzby1DXjBDcpOYN0E70eeAz47WXv30CjZQBwB/C6EMK69vdjjN8EqiGEbQPsXxnIW9fM9umN3HDljIFAGlDfLYMY4/8ACCF8dNk/vQSYb26zGEL4B6Da/n7TPLAVONTN/qamNgBQrU72e8iFkPfyveaSaf5s9yyLizUmJsZ4zSXTXR9zp+32zx5l38Gn2XnReeyY2Zz04Q5N3s/bICxbMbXqzG6tGQxCCG8GPrns7f0xxutX+Eilw+sajVZIvcP7XVlYOMHU1AaOHHmm248UTrU6mfvyTZ2zjvf9wumumalz1nV1zJ3Klqcup0EU4bz1y7IVU7U6earO7NaawSDGeBdwVw/HMQdcABwOIUwAk8ACcBjYAhxsbncB8GQP36ucSGrtoU5dTkUMBlIZpDG09M+BdzT/fiONZPLJ9vdDCFcDz8cYu+oiUjm53LWUH2lMOvsw8LkQwneBvwfe2nz/d4FPN99/AXh7CvtWgTgaSMqPgYNBa35B2+ujwD/vsN3zwC2D7k/lMszlrvMyP0LKI5ejUG6kWVmXJVktpcVgoFxIu7I2WS2tzrWJlAtpT2YzWS2tzpaBciHtdYZMVkurMxhoYEn09fdTWfe6X5/NLK3MYKCBJNnX30tlbUJYSpY5Aw0kq4Xr8rZgnlR0BgMNJKvErAlhKVl2E2kgWSVmTQhLyTIYaGDd9vUnPanMhLCUHIOBhmL/7FETvlKOmTPQUOw7+LQJXynHDAYaip0XnWfCV8oxu4k0FDtmNpvwlXLMloESdWDuOLv2zHJg7vgZ7++fPdpVIFjp8yu9LykZtgyUmJVmBR+YO84n7nyIk4urJ49X+7zJZyldtgyUmJVmBcdDx1hcXDt5vOrnTT5LqTIYKDErzQoO2zYxMbF28njVz5t8llJVqdfrWR/DWmaAJxYWTjA1tYEjR57J+nhSU61OFr587RPLgFN/33Tueh7cO9dVzqBTbiHPj6wsw3lbiWUrpmp1kladCbwcmF3rM+YMlKjWrODl/fwf+1dXccOVM11/vtP7wKkuorwFBKnoDAZKxfJ+/n0Hn+a6S7b0/X0mkaV0mTNQKpb38++86LyBvs8kspQuWwZKxfJVRXfMbB6ofzbtx2JKo85goNQkuaqoS1ZL6TIYKDM+w1jKD4OBMmFCWMoXE8jKhAlhKV8MBsqEs4qlfLGbSJkwISzli8FAmTEhLOWH3USSJIOBJMlgIEnCYCBJwmAgSaIYo4nGAcbGKrT/WVZlLp9lKybLVkxtZRvvZvsiPOnsauD+rA9CkgrqGuBba21UhGDwIuByYB5YyvhYJKkoxoEtwF8BL6y1cRGCgSQpZSaQJUkGA0mSwUCShMFAkoTBQJKEwUCShMFAkkQxlqMAIIRwM/AhYB3wOzHGT2V8SAMLIfwIsBv4pzHG2RDC9cDtwIuBL8YYP5TpAfYphPAR4C3Nl7tijB8oUdluBd4E1IHPxhhvL0vZWkIInwDOizG+M4RwKfAZ4EeAbwLvjjEuZnqAfQoh3AucD5xsvvXLwEWUoF4JIfwz4CPAOcDdMcZf7fW6LETLIIQwDdxGY2mKS4F3hRB+LNujGkwI4QoaU8Qvbr5+MfBHwOuBHwUuDyG8Lrsj7E/zAnwt8Eoa5+qyEMJNlKNs/wj4aeAS4FXAe0MIP0EJytYSQvgZ4Ja2tz4PvCfGeDFQAX4pkwMbUAihQuO39hMxxktjjJcChylBvRJCuBD4A+DnaVybP9m8Bnu6LgsRDIDrga/HGI/GGJ8FvkTj7qzIfgn4FeDJ5utXA4/FGJ9o3nl9HnhzVgc3gHng38cYfxBjPAn8DY0fYeHLFmP8BvCPm2U4n0bL+lxKUDaAEMJmGpXjx5qvXwa8OMb4YHOTz1HQsgGh+efdIYS/DiG8h/LUK2+gced/uPmbuxF4jh6vy6J0E72ERiXTMk+j8iysGOMvAoTQukY7lnHrkA9rYDHG77b+HkJ4BY3uot+lBGUDiDGeDCH8OvA+4C5Kct6aPg18EHhp83WZyrYJ+BrwXhpdQvcBX6Qc9cp24AchhD8FtgFfAb5Lj+euKC2DMRp9tC0VoJbRsaSlVGUMIfw48JfA+4HHKVHZYowfAao0Ks2LKUHZQgi/CPxdjPFrbW+X5pqMMe6JMb4jxng8xvg08FngVspRvgkarZx/CVwJXAFcSI9lK0rL4DCNZVhbLuB090pZHKaxwmBLYcsYQrgK+F/Ar8UY72z2tRe+bCGEHcDZMcaHY4zPhRC+TKNboX013UKWjUbXwpYQwsPAZmADjcqk8OcNIIRwNfCitmBXAWYpR/meAu6JMR4BCCH8MY0uoZ6uy6IEg3uAj4YQqsCzwBuBd2V7SIn7NhBCCNuBJ4CbaSSACiWE8FLgT4AbY4xfb75dirLRuNv69WbFUqeRnPs08PGily3G+E9afw8hvBO4Lsb4L0IIj4QQrooxPgC8HfhqVsc4oHOBW0MIP0Wjm+gW4G3A50tQr3wF+O8hhHOBZ4DX0ch//MderstCdBPFGOdo9GXeCzwM/M8Y43eyPapkxRifB95J4476UWA/jRNaNO8DzgZuDyE83LzTfCclKFuM8c+BXcBDwP8BdscY76QEZVvFW4FPhhD202gt/NeMj6cvMcavcOa5+6NmgCt8vRJj/DbwWzRGJz4K/C3w+/R4Xfo8A0lSMVoGkqR0GQwkSQYDSZLBQJKEwUCShMFAkoTBQJKEwUCSBPx/Y2rVVBuKKowAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "data = np.loadtxt(\"motor.dat\")\n",
    "X = data[:, 0][:, None]\n",
    "Y = data[:, 1]\n",
    "plt.plot(X, Y, '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution:**\n",
    "<br></br><br></br><br></br><br></br><br></br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part B - Hierarchical Priors \n",
    "\n",
    "Specify priors on the model hyperparameters and estimate the full joint posterior over the model weights and model hyperparameters. At a minimum, specify appropriate priors for the  prior precision of the weights vector, and the likelihood noise. Use the estimated posterior to get the posterior predictive distribution over test inputs. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution:**\n",
    "<br></br><br></br><br></br><br></br><br></br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part C - Heteroscedastic regression\n",
    "\n",
    "So far, throughout this course, you have seen likelihood models that assume that the noise level $\\sigma$ is independent of the input $x$. This is known as *homoscedasticity* - the assumption that errors in a regression model are indepedent of the inputs. Consider the following likelihood model with input dependent noise:\n",
    "$$\n",
    "y \\sim \\mathcal{N} ( y | \\mathbf{w}^T \\phi(x) , \\sigma(x)^2),\n",
    "$$\n",
    "where the likelihood noise depends on the input. \n",
    "Approximate $\\log \\sigma$ as a GLM of your choice - $\\log \\sigma = \\phi_{\\sigma}(x)^T \\mathbf{w_{\\sigma}}$ to model the dependence of the likelihood noise to the input. \n",
    "Develop the `PyMC3` model to express the heteroscedastic model and estimate the joint posterior over all parameters and hyperparameters. \n",
    "The parameters that you need to infer will include the weights of the output GLM $\\mathbf{w}$, the weights of the noise GLM model $\\mathbf{w}_{\\sigma}$, the precision over $\\mathbf{w}$ and the precision over $\\mathbf{w}_{\\sigma}$ and any additional hyperparameters you might have in your model.\n",
    "Use the estimated posterior to show the posterior predictive distribution over test inputs. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution:**\n",
    "<br></br><br></br><br></br><br></br><br></br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
